Pontos de corte em mapas

Ao representar uma variável em mapa (exemplo: número de casos, ou incidência de casos de uma doença), nem sempre representá-la de forma contínua pode ser informativo. Em distribuções desbalanceadas e/ou que possuem valores aberrantes (valores próximos de zero e valores muito mais altos que os demais), a escala de cores para o mapa pode acabar chamando atenção somente para valores extremos.

Vamos tomar como exemplo a incidência de casos prováveis de dengue por município no estado de São Paulo no ano de 2015. A partir do número de casos notificados no SINAN e da população estimada para o ano em questão, podemos calcular a incidência e representá-la, de forma contínua, em um mapa como na Figura 18.

Figura 18: Incidência de casos prováveis de dengue por município no estado de São Paulo, 2015.

No mapa, nota-se alguns municípios que chamam mais atenção pela cor mais escura - ou seja, maior incidência. Alguns outros pontos de concentração da doença são marcados no estado, mas, no geral, a maioria dos municípios é pintada em tom mais claro.

A visualização da variável em sua forma contínua acaba sendo útil em poucos casos, principalmente em situações ligadas à vigilância em saúde. Tradicionalmente, trabalhamos com dados de contagem de casos - e, na maioria dos cenários, a distribuição dessa variável é desbalanceada. Poucos municípios possuem muitos casos enquanto muitos municípios possuem poucos casos. Além disso, municípios com populações pequenas acabam tendo a incidência inflada por conta de um denominador pequeno. Isso faz com que poucos municípios se destaquem dos demais na escala de cores utilizando a variável de forma contínua, e torna difícil a distinção do cenário entre o restante.

Dessa forma, torna-se importante a separação dessa variável em faixas - através da definição de pontos de corte. Contudo, a escolha da estratégia utilizada para definição desses pontos importa e pode impactar e muito no resultado da visualização obtida.

Definindo pontos de corte

Observe os mapas abaixo na Figura 19. Eles representam a mesma informação - incidência de casos prováveis de dengue por município no estado de São Paulo, em 2015. O que há, portanto, de diferente entre eles?

Figura 19: Formas de representação da incidência de casos prováveis de dengue no estado de São Paulo, 2015.

Podemos observar que, a depender da estratégia de definição dos pontos de corte utilizada, a interpretação que temos de cada mapa pode ser diferente. Em alguns mapas (A e C), a situação no estado parece pior, com mais municípios pintados em tons escuros, levando a pensar em uma maior incidência. No mapa B, no entanto, a impressão é exatamente a oposta - a maioria dos municípios está em tons claros, o que nos leva a interpretar que o estado em geral está com baixa incidência da doença. Vamos entender cada uma dessas situações.

Quebras arbitrárias

O primeiro caso, representado na Figura 19A, não utiliza nenhum método específico para estabelecimento de pontos de corte da incidência. Os pontos são definidos arbitrariamente, e podem ou não ser baseados em uma referência existente sobre o tema. Em um departamento de vigilância, pode-se ter estabelecido que uma incidência de 50 casos por 10.000 habitantes é um patamar controlado, e acima de 100 casos por 100.000 tem-se um cenário de preocupação. Assim, o analista pode-se basear nesses pontos pré-definidos na rotina de vigilância para gerar seus mapas. Contudo, não é garantia de que a visualização seja adequada aos dados: pode ser que muitos municípios se concentrem em uma mesma faixa, por exemplo, trazendo dificuldades na distinção entre esses municípios (Figura 20).

Figura 20: Distribuição dos municípios em cada faixa de incidência, definida por quebras arbitrárias.

Quebras regulares

A segunda estratégia de quebras utilizada é a das quebras regulares. Nela, considera-se o intervalo entre o mínimo e o máximo da distribuição, e divide-se esse intervalo em partes iguais. Dessa forma, cada faixa das quebras possui a mesma amplitude.

Repare, na Figura 19B, que a distância entre um ponto de corte e outro é de aproximadamente 296 para todos os intervalos.

Essa estratégia facilita a comunicação, mas é particularmente útil quando a distribuição dos dados se aproxima de uma distribuição uniforme. No cenário da vigilância em saúde, esse é raramente o caso - há muitos valores concentrados em uma faixa da distribuição e poucos valores distribuídos entre outras. Isso causa o efeito de que um número excessivo de municípios pode cair na mesma categoria, enquanto outras categorias podem possuir um ou dois municípios (Figura 21). Pode ser, até, que faixas de incidência não possuam nenhum município dentro dela.

Figura 21: Distribuição dos municípios em cada faixa de incidência, definida por quebras regulares.

Vê-se agora que a distância entre a mudança de cores é sempre a mesma, o que é agradável visualmente. Contudo, a vasta maioria dos municípios se concentra dentro de uma mesma cor. Dessa forma, a estratégia falha em distinguir esses municípios com incidência mais baixa entre eles.

Quebras quantílicas

Outra estratégia consiste em determinar pontos de quebra a partir de quantis da distribuição. As quebras quantílicas tentam dividir os municípios em grupos de forma com que cada grupo contenha aproximadamente a mesma quantidade de municípios. Por exemplo, ao dividir em 5 grupos, teríamos cada grupo com aproximadamente 20% do total de observações.

Figura 22: Distribuição dos municípios em cada faixa de incidência, definida por quebras quantílicas

Dessa forma, “forçamos” os grupos a conterem aproximadamente a mesma quantidade de municípios. Essa técnica garante que os mapas terão representação visível de todos os grupos; contudo, vemos que mesmo municípios com incidências baixas (Figura 22, primeira barra) já estão sendo categorizados para o segundo grupo, enquanto o ponto de corte para o grupo de incidência mais alta fica muito mais baixo. Isso pode levar a falsas impressões nos mapas, já que fazemos com que necessariamente uma parcela igual esteja contida nos níveis mais altos e mais baixos de incidência.

Quebras naturais, ou quebras de Jenks

O método de classificação por quebras naturais (ou quebras de Jenks) parte da ideia de obter pontos de corte que tentam, ao mãximo, minimizar a variação intra-grupo e maximizar a variância entre os grupos. Ou seja, busca-se dividir o conjunto de dados em grupos semelhantes entre si e diferentes dos demais.

Figura 23: Distribuição dos municípios em cada faixa de incidência, definida por quebras naturais de Jenks.

Para este exemplo, as quebras naturais classificaram os municípios de forma equilibrada em relação às demais (Figura 23) - não houve uma hiper concentração dos municípios em uma só classe (como as quebras regulares), nem a definição de quebras que abrangiam muitos municípios diferentes dentro de si (como nas quebras quantílicas). Percebe-se que o primeiro grupo consta majoritariamente de municípios com baixa incidência (primeira barra no histograma), e o grupo mais alto inclui poucos municípios, mas que realmente apresentaram valores de incidência mais extremos na distribuição.

Comparando mapas no tempo

Outro ponto importante ao definir escalas e pontos de corte em mapas é a padronização dessas escalas ao comparar mapas em diferentes pontos no tempo.

Por exemplo, suponha que queremos comparar a incidência de dengue nos municípios do estado de São Paulo no período 2014-2017, utilizando as quebras quantílicas (Figura 24).

Figura 24: Incidência de dengue nos municípios São Paulo, ao longo de 2014-2017.

Os mapas gerados nos levam a interpretações sobre as áreas onde a doença se concentrou em cada ano. Podemos dizer que há certa similaridade entre os padrões espaciais, onde certas regiões do Norte e Noroeste paulista aparecendo repetidamente nas faixas de maior incidência. A impressão geral que temos é que, a distribuição da doença se deu de forma parecida ao longo dos anos. Contudo, vemos que os pontos de corte em cada mapa são criticamente diferentes: enquanto o tom mais escuro no mapa relativo à 2017 representa uma incidência maior que 5,67, o mesmo tom no mapa de 2015 representa incidência superior a 393,97.

Assim, ao comparar mapas no tempo, faz-se necessária a utilização de uma escala única para permitir a comparação das taxas não somente no espaço, mas entre os anos. Veja:

Figura 25: Incidência de dengue nos municípios São Paulo, ao longo de 2014-2017, utilizando escalas de cor independentes a cada ano.

Agora, utilizamos a mesma técnica de classificação a partir de quebras quantílicas, mas para todo o período de estudo para determinação das quebras. Isso permite a comparação dos cenários ao longo do tempo: além de identificarmos certos padrões espaciais, vemos que o ano de 2015 foi muito mais crítico em termos de incidência de dengue do que os demais anos - principalmente 2017, quando a maioria dos municípios ficou no primeiro grupo da escala de cores.

Tamanho heterogêneo das áreas

Como mencionado, as delimitações no território frequentemente são derivadas de divisões políticas e administrativas. Entre outros efeitos, isso pode levar à geração de áreas pequenas em extensão territorial, porém com alta densidade populacional; enquanto áreas extensas concentram pouca população. Isso pode afetar nossas visualizações: as áreas maiores, além de chamarem maior atenção aos nossos olhos, estão sujeitas a uma inflação das taxas devido ao baixo denominador (população). Vamos conferir, por exemplo, a análise da taxa de alfabetização dos municípios brasileiros, segundo o Censo 2022 (Figura 26).

As áreas de principal destaque no mapa são áreas consideravelmente grandes, e que apresentam piores indicadores.

Figura 26: Taxa de Alfabetização (%) nos municípios brasileiros, segundo o Censo 2022.

Vamos comparar esse padrão espacial com a densidade populacional nesses municípios (Figura 27).

Figura 27: Densidade populacional (habitantes/km²) nos municípios brasileiros, segundo o Censo 2022.

Vemos que além de grandes, as áreas de destaque no primeiro mapa possuem baixa densidade populacional. Como a população é frequentemente o denominador nos indicadores que analisamos, quando a população é pequena essas taxas tendem a ser instáveis. Portanto, municípios de população pequena (e que, como vimos, podem ser grandes em termos de extensão territorial) estão sujeitos a uma elevação artificial nas taxas devido a essa baixa população, e acabam mascarando demais padrões espaciais quando olhamos o mapa completo.

Problema da área modificável (MAUP)

Outro problema comum é que as taxas e indicadores dependem da forma com que as áreas são delimitadas e agregadas. E se, de repente, essas delimitações são alteradas, resultados diferentes podem ser obtidos. Esse problema é conhecido como Problema da Área Modificável (Modifiable Areal Unit Problem - MAUP).

Vamos tentar ilustrar esse efeito a partir de dados hipotéticos de casos de uma doença sobre um conjunto de setores censitários (Figura 28).

Figura 28: Número de casos de uma doença hipotética sobre um conjunto de setores censitários.

Agora, vamos supor que esses dados são agregados para uma unidade espacial maior (como bairro) para serem disponibilizados. A forma com que essas delimitações de bairro são definidas a partir dos setores censitários pode mudar a visualização da incidência da doença que obtemos. Vamos comparar duas delimitações distintas na Figura 29:

Figura 29: Comparação da incidência da doença hipotética em duas divisões de território diferentes (A e B) a partir dos setores censitários.

Vemos que as duas divisões, embora passem a mesma informação de distribuição da doença, mostram impressões ligeiramente diferentes de onde ela se concentra. Além disso, taxas de incidência consideravelmente mais altas são obtidas a depender do zoneamento. Como o processo de agregação pode se dar em vários níveis a partir dos dados referenciados, deve-se ter atenção para esses possíveis efeitos de atenuação ou de amplificação das taxas de incidência a depender dos níveis e divisões no processo de agregação.

Gerrymandering

A variação de proporções e taxas de acordo com a mudança dos limites territoriais pode ser usado intencionalmente e com a finalidade de alterar resultados importantes - como é feito no contexto eleitoral estadunidense. Em tal país, parte do processo eleitoral envolve contar qual partido é ganhador em cada distrito (voto distrital) e a partir disso, contar quantos distritos cada partido “tem”. A forma com que os distritos são delimitados é extremamente importante, pois pode alterar os resultados obtidos no final das contas.

Observe o exemplo a seguir na Figura 30, onde temos 50 pessoas e 60% delas votaram no partido laranja, e os 40% restantes no partido roxo.

Figura 30: Exemplo de Gerrymandering - A depender da divisão dos distritos, o resultado da contagem é alterado.

Perceba que, dependendo da forma com que são traçados os distritos, temos uma mudança significativa na contagem de distritos que cada partido possui: no terceiro cenário, temos todos os distritos com vitória do partido laranja. Por outro lado, mesmo com menor proporção total de votos no total da população, a quarta configuração faz com que o partido roxo tenha mais distritos - 3 contra 2 do laranja. Esse procedimento de calcular as divisões distritais com a finalidade de alterar os resultados obtidos no processo eleitoral ocorre nos Estados Unidos, e é chamado de Gerrymandering.

Malhas geográficas com geobr

Agora vamos visualizar nossos dados no espaço!

Tradicionalmente, o carregamento de arquivos espaciais para dentro do ambiente R é feito através da leitura de shapefiles (arquivos que incorporam a dimensão espacial aos dados), e usualmente através da função st_read(). Essa função pode ser utilizada da seguinte forma:

malha_exemplo <- st_read("caminho/para/o/shape/malha_br.shp").

Contudo, com a ampla comunidade colaborativa do R, certas operações têm se tornado mais práticas. Hoje, o pacote geobr, desenvolvido pelo IPEA, reune malhas territoriais oficiais de diferentes fontes, como IBGE, CEMADEN e INEP em um repositório único que pode ser acessado diretamente via R, sem necessidade de baixar nenhum shapefile manualmente. Vamos ver como isso funciona:

# se não estiver instalado, rodar:
install.packages("geobr")
library(geobr)

O pacote possui funções específicas para cada nível de agregação espacial, todas começando com read_*. Por exemplo, se quiséssemos uma malha dos estados brasileiros, poderíamos usar read_state(). Para bairros, read_neighborhood(). No caso, iremos recuperar a malha de municípios de São Paulo, e, portanto, utilizaremos read_municipality():

# carregando a malha do estado de São Paulo
malha_sp <- read_municipality("SP", showProgress = F)
head(malha_sp)

#> Simple feature collection with 6 features and 4 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -51.17838 ymin: -23.03648 xmax: -46.54881 ymax: -21.19701
#> Geodetic CRS:  SIRGAS 2000
#>   code_muni              name_muni code_state abbrev_state
#> 1   3500105             Adamantina         35           SP
#> 2   3500204                 Adolfo         35           SP
#> 3   3500303                  Aguaí         35           SP
#> 4   3500402         Águas Da Prata         35           SP
#> 5   3500501       Águas De Lindóia         35           SP
#> 6   3500550 Águas De Santa Bárbara         35           SP
#>                             geom
#> 1 MULTIPOLYGON (((-51.09093 -...
#> 2 MULTIPOLYGON (((-49.69668 -...
#> 3 MULTIPOLYGON (((-47.01254 -...
#> 4 MULTIPOLYGON (((-46.73069 -...
#> 5 MULTIPOLYGON (((-46.635 -22...
#> 6 MULTIPOLYGON (((-49.28903 -...

Vemos que os dados estão dispostos de forma similar a uma base de dados tradicional, possuindo agora uma coluna espacial: a coluna geom. Ela contém as informações dos polígonos referentes a cada um dos municípios. É importante notar também o Sistema de Coordenadas (CRS) listado - nesse caso, o SIRGAS 2000.

Podemos verificar rapidamente o contorno dos municípios:

plot(malha_sp$geom)

Vemos que em ambos os objetos (dados de dengue e malha espacial) possuímos a coluna referente ao código do município, e, portanto, a usaremos para juntar as duas informações - as informações epidemiológicas e espaciais.

Podemos realizar essa operação através do comando left_join(). Perceba que o código do município é chamado de code_muni na malha espacial e cod_ibge no conjunto de dados de dengue:

malha_sp_dengue <- malha_sp %>% 
  left_join(
    dengue_2015 %>% select(cod_ibge, casos, pop, inc),
    by = c("code_muni" = "cod_ibge")
  )

head(malha_sp_dengue)

#> Simple feature collection with 6 features and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -51.17838 ymin: -23.03648 xmax: -46.54881 ymax: -21.19701
#> Geodetic CRS:  SIRGAS 2000
#>   code_muni              name_muni code_state abbrev_state casos   pop
#> 1   3500105             Adamantina         35           SP  1140 34285
#> 2   3500204                 Adolfo         35           SP   121  3903
#> 3   3500303                  Aguaí         35           SP  2487 32192
#> 4   3500402         Águas Da Prata         35           SP   336  7558
#> 5   3500501       Águas De Lindóia         35           SP   165 17690
#> 6   3500550 Águas De Santa Bárbara         35           SP    16  6313
#>         inc                           geom
#> 1 332.50693 MULTIPOLYGON (((-51.09093 -...
#> 2 310.01793 MULTIPOLYGON (((-49.69668 -...
#> 3 772.55219 MULTIPOLYGON (((-47.01254 -...
#> 4 444.56205 MULTIPOLYGON (((-46.73069 -...
#> 5  93.27304 MULTIPOLYGON (((-46.635 -22...
#> 6  25.34453 MULTIPOLYGON (((-49.28903 -...

Agora sim! Temos as informações espaciais (coluna geom) e epidemiológicas juntas em um objeto. A partir disso, podemos construir nosso primeiro mapa de incidência de dengue em São Paulo no ano de 2015. Vamos utilizar a função ggplot() que possui uma função específica para visualização de mapas: a geom_sf(). Precisamos especificar duas estéticas (aes()): a geometry, que recebe o nome da coluna que possui as informações espaciais; e fill, que recebe o nome da variável pela qual o mapa será colorido:

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc), linewidth = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, 2015", fill = "Incidência (por 10.000 hab.)") + 
  theme_void()

Podemos alterar a paleta de cores através dos comandos scale_fill_*, como scale_fill_distiller() e scale_fill_viridis_c() para variáveis contínuas e scale_fill_brewer() e scale_fill_viridis_d() para variáveis discretas ou categorizadas.

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc), linewidth = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, 2015", fill = "Incidência (por 10.000 hab.)") + 
  scale_fill_distiller(palette = "Reds", direction=1) +
  theme_void()

Pronto! Temos nosso primeiro mapa de incidência de dengue. Conforme discutido anteriormente, o mapa representa uma taxa contínua, que nem sempre é a melhor opção para visualizar os padrões espaciais - ainda mais quando a distribuição da variável não é simétrica.

Vamos definir, portanto, as quebras - pontos de corte para divisão da incidência no mapa, conforme estratégias que vimos anteriormente.

# quebras pré-definidas (arbitrárias)
quebras_arbitrarias <- c(0, 10, 50, 100, 250, 500, 2000)

malha_sp_dengue <- malha_sp_dengue %>% 
  mutate(
    inc_cat = cut(inc, breaks = quebras_arbitrarias, include.lowest=T)
  )

E, ao realizar o gráfico, usamos inc_cat como variável e scale_fill_brewer() para especificar as cores:

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc_cat), linewidth = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, 2015", 
 fill = "Incidência (por 10.000 hab.)") + 
  scale_fill_brewer(palette = "Reds", direction=1) +
  theme_void()

Vamos explorar as outras formas de definir as quebras (fique à vontade para trocá-las no mapa e comparar as diferenças).

# quebras regulares
quebras_regulares <- seq(0,
                         max(dengue_2015$inc), 
                         length.out=7) # length.out define o número de categorias + 1

# quebras quantilicas
quebras_quantilicas <- quantile(
                         dengue_2015$inc, 
                         probs = seq(0, 1, length.out=7))

Para as quebras de Jenks (ou quebras naturais), é necessário ter instalada uma implementação do algoritmo de Jenks. Uma opção é através da biblioteca BAMMtools:

install.packages("BAMMtools")
library(BAMMtools)

# quebras naturais ou quebras de Jenks

quebras_jenks <- BAMMtools::getJenksBreaks(
  dengue_sp$inc, k=7
)

quebras_regulares

#> [1]    0.0000  295.5472  591.0944  886.6416 1182.1888 1477.7360 1773.2832

quebras_quantilicas

#>         0%  16.66667%  33.33333%        50%  66.66667%  83.33333%       100% 
#>    0.00000   33.81285   79.98215  135.09475  232.76623  393.97020 1773.28316

quebras_jenks

#> [1]    0.00000   59.01676  174.00383  329.79113  556.57492  933.26643 1773.28316

malha_sp_dengue <- malha_sp_dengue %>% 
  mutate(
    inc_cat2 = cut(inc, breaks = quebras_jenks, include.lowest=T)
  )

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc_cat2) , linewidth = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, 2015", 
) +
  theme_void()

Pronto! Através da junção de dados epidemiológicos com malhas espaciais construímos nossas primeiras visualizações em mapas com dados de área, alternando paletas e maneiras de definir os pontos de corte.

Vamos voltar agora para o conteúdo teórico do curso!

Efeitos de primeira ordem - tendência

Quando a dependência espacial sobre os dados se dá devida a efeitos de uma variação de larga escala (ou estrutural), ela é caracterizada como um efeito de primeira ordem. Esses efeitos são observados quando, por exemplo, há tendência de concentração da doença em regiões litorâneas: há uma determinação geográfica subjacente que induz essa tendência, levando a maiores valores de incidência em uma determinada região. Outro caso pode acontecer quando há uma tendência clara de aumento de incidência de uma doença no sentido Norte - Sul.

Vejamos um exemplo de efeito de primeira ordem. Retornando ao nosso conjunto de dados hipotético, vamos supor que uma doença X tenha a seguinte distribuição em um município (suponha população aproximadamente igual entre as áreas):

Figura 31: Exemplo de efeito espacial de primeira ordem: conforme vamos em direção ao leste, os casos da doença X parecem aumentar.

Vemos na Figura 31 que a distribuição da doença se dá de forma mais intensa conforme andamos em sentido Leste (ou Sudeste). Pode-se considerar que há, portanto, uma tendência “macro” no processo, em que o mesmo aumenta de intensidade ao avançar em determinado sentido, caracterizando um efeito de primeira ordem. Geralmente, estes efeitos se dão devido a alguma característica ambiental, como maior presença de áreas verdes, maior proximidade com corpos fluviais, maior/menor temperatura, umidade, pressão, entre outras.

Efeitos de segunda ordem - dependência local

Por outro lado, efeitos de segunda ordem (ou efeitos meso, ou efeitos locais) caracterizam uma dependência de uma região sobre seus vizinhos mais próximos. Um exemplo clássico é quando observamos “clusters” de uma doença infecciosa em algumas regiões: se há uma alta de casos em um determinado município, é muito provável que haja também uma alta de casos nos municípios vizinhos, seja pela própria dinâmica infecciosa da doença ou por características ambientais locais propícias para o desenvolvimento de vetores.

Vamos agora nos voltar ao exemplo dos casos de dengue em São Paulo, e verificar a incidência da distribuição de casos em 2020:

dengue_2020 <- dengue_sp %>% 
  filter(ano == 2020)

malha_sp_dengue <- malha_sp %>% 
  left_join(
    dengue_2020 %>% select(cod_ibge, casos, pop, inc),
    by = c("code_muni" = "cod_ibge")
  )

Vamos optar pela classificação das faixas pelas quebras de Jenks:

quebras_jenks <- BAMMtools::getJenksBreaks(
  dengue_2020$inc, k=7
)

malha_sp_dengue <- malha_sp_dengue %>% 
  mutate(
    inc_cat = cut(inc, breaks = quebras_jenks, include.lowest=T)
  )

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc_cat), linewdith = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, em 2020", 
       subtitle = "Quebras de Jenks", 
       fill = "Incidência (por 10.000 hab.)") + 
  scale_fill_brewer(palette = "Oranges", direction=1) +
  theme_void()

Figura 32: Incidência de dengue nos municípios de São Paulo, 2020.

Vemos na Figura 32 que há certos “grupos de concentração” da doença no estado. Pode-se listar ao menos quatro regiões em que há um foco de incidência maior entre municípios próximos. Não vemos, no geral, uma tendência clara ao longo de todo o estado como no exemplo anterior. Nesse caso, temos aglomerados da doença nas regiões Oeste e Noroeste do estado, assim como no Sul e Leste também. Esse processo de dependência espacial se aproxima, portanto, de um processo de segunda ordem.

Medidas de proximidade em dados de área

Vimos que as definições de autocorrelação espacial relacionam a variável de interesse entre as áreas próximas ou vizinhas. Portanto, para identificar, testar e quantificar a presença de autocorrelação nos dados de interesse, temos que definir quais áreas são próximas das outras a partir de medidas de proximidade. Essas podem se basear exclusivamente nas fronteiras territoriais, ou em outros fatores como distância (exemplo: distância entre os centroides, ou distância entre suas sedes) ou fluxo (exemplo: fluxo pendular) entre os municípios. Essa escolha se dará a partir dos mecanismos envolvidos no estudo do agravo em questão; há bastante sentido, por exemplo, em considerar o fluxo de deslocamento entre municípios para relacionar espacialmente a distribuição de uma doença como a covid-19. Dessa forma, municípios que tem maior fluxo de pessoas entre si serão mais “próximos” do que municípios geograficamente próximos, mas que não possuem tanto fluxo. Contudo, a disponibilidade de dados sobre esses fatores pode ser um fator limitador, o que nos leva a optar por abordagens mais simples. Vejamos a Figura 33. Neste caso, calculamos a proximidade com base na distância entre os centroides das áreas. Quanto mais perto duas áreas são, mais forte (pois mais próxima) é sua ligação (representada pelas linhas azuis).

Figura 33: Exemplo de cálculo de proximidade entre as áreas do exemplo hipotético.

Essas distâncias serão incorporadas a um elemento essencial ao tratar de autocorrelação espacial: a matriz de vizinhança, que veremos a seguir.

Matriz de vizinhança

Uma matriz de vizinhança é uma forma que temos de relacionar as áreas entre si, indicando quais possuem conexões (vizinhança) com as outras e quais os pesos dessas conexões. Geralmente temos uma matriz $n \times n$, onde temos todos os municípios nas linhas e os mesmos municípios nas colunas. Cada célula então, de uma linha $l$ com uma coluna $c$ indica se o município $l$ e $c$ tem relação de vizinhança - e qual o peso dessa relação. Os pesos podem ser fixos (como 0 para uma relação de não vizinhança e 1 para uma relação de vizinhança) ou dependerem de algum outro fator, como a distância entre os municípios.

Vamos considerar o seguinte exemplo, de um recorte de alguns municípios do estado de São Paulo (Figura 34).

Figura 34: Recorte de 7 municípios do estado de São Paulo.

Uma possível representação das relações entre os municípios exibidos e sua vizinhança seria:

Vemos que a primeira linha, que representa o município de Arujá, possui valor 1 nas colunas Itaquaquecetuba e Mogi das Cruzes, pois são os municípios com quem faz fronteira. Mogi das Cruzes, por sua vez, possui valor 1 para quase todos os municípios listados, com exceção de Salesópolis. A matriz exibida é chamada de matriz de vizinhança por contiguidade binária; pois define se um município é vizinho de outro se compartilham fronteira, e caso exista, o peso é atribuído é 1. Caso contrário, o peso é 0. Este tipo de matriz de vizinhança é chamado de vizinhança por contiguidade.

As estratégias mais comuns para definição da matriz de vizinhança são:

• Por contiguidade: Como no exemplo acima, a célula $l,c$ assumirá valor 1 se o município $l$ e o município $c$ compartilharem fronteira e 0 caso contrário;

• Por $k$ vizinhos mais próximos: a célula $l,c$ será 1 se $c$ estiver entre os $k$ municípios mais próximos de $l$. Pode-se considerar, por exemplo, a distância entre as sedes ou centroides dos municípios;

• Por distância: a célula $l,c$ será 1 se $c$ estiver até a uma distância $d$ pré-definida de $l$. Caso contrário, será 0.

Vamos implementar cada uma dessas estratégias em nosso mapa do estado de São Paulo.

Por contiguidade

O pacote spdep traz diversas operações envolvendo correlação e dependência espacial implementadas, e, portanto, o utilizaremos para definir nossas matrizes de vizinhança. Caso não esteja instalado, precisaremos rodar install.packages("spdep").

# se não estiver instalado, rodar:
install.packages("spdep")
library(spdep)

A library(spdep) no R é usada para análise de estatística espacial. Ela fornece funções para criar estruturas de vizinhança, calcular pesos espaciais e aplicar testes e modelos que mensuram a autocorrelação espacial.

A partir da malha de municípios de São Paulo, vamos definir nossa matriz de vizinhança. O pacote spdep possui diferentes funções de transformação nesse sentido, onde podemos partir de objetos espaciais como o que temos (malha_sp) e transformá-lo em objetos como de relação entre vizinhos. Uma dessas funções é a poly2nb (polígono (poly) para (to - 2) vizinhança (neighbours - nb)). Ela automaticamente define uma matriz de contiguidade a partir da presença de fronteiras na malha espacial em questão. Vamos testá-la e criar o objeto viz_sp, que utilizaremos mais adiante.

viz_sp <- poly2nb(malha_sp)
viz_sp

#> Neighbour list object:
#> Number of regions: 645 
#> Number of nonzero links: 3530 
#> Percentage nonzero weights: 0.8485067 
#> Average number of links: 5.472868 
#> 1 region with no links:
#> 233
#> 2 disjoint connected subgraphs

A função poly2nb(malha_sp) cria uma lista de vizinhança com base no polígono malha_sp.

• O objeto malha_sp deve ser um objeto sf ou SpatialPolygons.

• A função poly2nb() é do pacote spdep.

• Ela define vizinhos com base no compartilhamento de bordas (isto é, dois polígonos são vizinhos se eles dividem um lado).

• O objeto viz_sp representa o objeto resultante, que é da classe nb (neighbour list). Cada elemento da lista contém os índices dos polígonos vizinhos de cada polígono de malha_sp.

• viz_sp impresso no console vai mostrar para cada polígono, quais são os seus vizinhos.

Ótimo! Temos, então, um objeto do tipo “lista de vizinhos” (Neighbour list). Ao chamar o objeto, temos um resumo dessa matriz gerada: são 645 municípios, 3.526 links (ligações), tendo cada região em média 5.47 links (ou seja, 5 vizinhos).

Somos informados também, que uma região (de id 233) não possui links, ou seja, não possui vizinhos. Isso geralmente é um problema quando falamos de análise de dependência espacial. Para testar nossas hipóteses de autocorrelação espacial e utilizar outras técnicas relacionadas, é ideal que todas as regiões tenham ao menos um vizinho. Pode acontecer que, em caso de linhas ou regiões desconexas, ao gerar a lista de vizinhos automaticamente a partir da malha, alguma dessas regiões fique desconectada das demais. Mas lidaremos com isso mais à frente.

Vamos tentar visualizar nossa vizinhança obtida:

# obtendo as coordenadas dos centroides da malha

centroides_sp <- st_centroid(malha_sp)
centroides_coordenadas <- st_coordinates(centroides_sp)

# a partir delas, transformamos nossas relações de vizinhança em linhas:

viz_linhas <- nb2lines(nb = viz_sp, 
                       coords = centroides_coordenadas) %>% 
  st_as_sf()
st_crs(viz_linhas) <- st_crs(malha_sp)

ggplot(malha_sp, aes(geometry=geom)) + 
  geom_sf(fill = "lightblue") + 
  geom_sf(data = viz_linhas, aes(geometry=geometry)) + 
  ggtitle("Matriz de vizinhança por contiguidade, estado de São Paulo.") +
  theme_void()

Figura 35: Matriz de vizinhança por contiguidade dos munícipios do estado de São Paulo.

Esses códigos criam e desenham as relações de vizinhança entre os municípios do estado de São Paulo, representando essas relações com linhas conectando os centroides (pontos centrais) dos polígonos da malha.

• centroides_sp <- st_centroid(malha_sp): calcula o ponto central (centroide) de cada área (polígono) da malha malha_sp.

• centroides_coordenadas <- st_coordinates(centroides_sp): Extrai as coordenadas desses centroides para trabalhar com elas depois.

• viz_linhas <- nb2lines(nb = viz_sp, coords = centroides_coordenadas) %>% st_as_sf(): Converte a lista de vizinhos viz_sp em linhas, usando as coordenadas dos centroides. O resultado é transformado em um objeto sf (simple features).

• st_crs(viz_linhas) <- st_crs(malha_sp): Define o sistema de referência de coordenadas (CRS) do objeto viz_linhas para ser o mesmo que o da malha original malha_sp.

• Logo em seguida com a função ggplot, o mapa é plotado primeiro com a malha dos polígonos (malha_sp) com cor azul clara. E depois com as linhas de vizinhança sobre a malha, ligando áreas vizinhas.

Vemos na Figura 35 que todos os municípios que tocam suas fronteiras estão conectados. Os municípios da “borda” do estado tendem a ter menos vizinhos (alguns têm apenas um) e agora conseguimos visualizar o município sem nenhum vizinho. Trata-se de Ilhabela, no Litoral Norte de São Paulo, que justamente por ser uma ilha, não possui fronteira física com nenhum outro município na malha considerada.

Vamos, então, voltar ao nosso objeto de vizinhança viz_sp.

viz_sp

#> Neighbour list object:
#> Number of regions: 645 
#> Number of nonzero links: 3530 
#> Percentage nonzero weights: 0.8485067 
#> Average number of links: 5.472868 
#> 1 region with no links:
#> 233
#> 2 disjoint connected subgraphs

Somos informados que o município que não possui vizinhos é o 233. Vamos verificar se é realmente Ilhabela:

malha_sp$name_muni[233]

#> [1] "Ilhabela"

Pegando o nome do município que está na posição 233 do objeto malha_sp.

Sim, é Ilhabela. Vamos, portanto, conectá-la manualmente a um outro município. Como a chegada a Ilhabela é feita através de balsa partindo de São Sebastião, faz sentindo conectá-la a esse município.

Vamos verificar qual o índice que corresponde a São Sebastião:

which(malha_sp$name_muni == "São Sebastião")

#> [1] 567

Procura a posição (o número do índice) onde o nome do município é “São Sebastião” dentro da coluna name_muni do objeto malha_sp.

Agora que sabemos o índice dos dois municípios, podemos modificar nosso objeto viz_sp para conectá-los:

Incluindo São Sebastião como vizinho de Ilhabela:

viz_sp[[233]] <- 567L

Está forçando a região número 233 (que antes não tinha vizinhos) a ter como vizinho o município de número 567.

Note que utilizamos a notação 567L ao invés de somente 567. Isso se deve porque a estrutura de vizinhança aceita somente índices inteiros. No ambiente R, quando trabalhamos com números sem informar explicitamente seu tipo, o ambiente assuma que é um double (número real de precisão dupla):

class(567)

#> [1] "numeric"

class(567L)

#> [1] "integer"

is.double(567)

#> [1] TRUE

is.double(567L)

#> [1] FALSE

O comando viz_sp[[233]] <- 567L altera manualmente a lista de vizinhança espacial viz_sp, atribuindo ao município de índice 233 o número 567 como seu vizinho. A letra L após o número indica que ele é tratado como um inteiro (integer) no R, o que é necessário porque listas de vizinhança esperam índices inteiros. Essa modificação é útil para corrigir problemas de regiões isoladas, garantindo que todas as áreas tenham pelo menos um vizinho para análises espaciais posteriores.

Logo, precisamos ter essa rigorosidade ao alterar os vizinhos senão o pacote lançará um erro. Agora, vamos fazer o inverso (que também é necessário!): incluir Ilhabela como vizinho de São Sebastião:

viz_sp[[567]] <- c(viz_sp[[567]], 233L)

O comando viz_sp[[567]] <- c(viz_sp[[567]], 233L) está atualizando a lista de vizinhança do município 567, adicionando o município 233 como seu novo vizinho. Ou seja, ele pega os vizinhos que 567 já tinha (viz_sp[[567]]) e junta (c(...)) o número 233 (especificado como inteiro com L), garantindo que a relação de vizinhança seja bidirecional: se 233 tem 567 como vizinho, então 567 também passa a ter 233 como vizinho.

E agora, se repetirmos o comando do gráfico realizado anteriormente, vemos a inclusão da ligação entre Ilhabela e São Sebastião (Figura 36):

viz_linhas <- nb2lines(nb = viz_sp, 
                       coords = centroides_coordenadas) %>% 
  st_as_sf()
st_crs(viz_linhas) <- st_crs(malha_sp)

ggplot(malha_sp, aes(geometry=geom)) + 
  geom_sf(fill = "lightblue") + 
  geom_sf(data = viz_linhas, aes(geometry=geometry)) + 
  labs(title="Matriz de vizinhança por contiguidade, estado de São Paulo.", 
       subtitle = "Após inclusão de Ilhabela") +
  theme_void()

Figura 36: Matriz de vizinhança por contiguidade modificada dos munícipios do estado de São Paulo.

O código acima transforma a lista de vizinhança viz_sp em linhas que conectam os centroides das regiões (centroides_coordenadas) usando a função nb2lines, depois converte essas linhas em um objeto espacial (sf) com st_as_sf(). Em seguida, ajusta o sistema de coordenadas das linhas para ser igual ao da malha malha_sp com st_crs(). Por fim, cria um mapa com ggplot2, desenhando as áreas (malha_sp) preenchidas de azul claro e sobrepondo as linhas de vizinhança (viz_linhas), adicionando título, subtítulo e removendo elementos gráficos extras (theme_void()).

Verificando o objeto de vizinhança novamente:

viz_sp

#> Neighbour list object:
#> Number of regions: 645 
#> Number of nonzero links: 3532 
#> Percentage nonzero weights: 0.8489874 
#> Average number of links: 5.475969 
#> 2 disjoint connected subgraphs

Agora sim! Todos os municípios estão conectados, e não temos nenhum município sem link.

Por vizinhos mais próximos

A segunda estratégia possível que vimos é a definição dos vizinhos a partir dos $k$ municípios mais próximos a um município específico. Repare essa estratégia no exemplo, considerando os dois municípios mais próximos ($k=2$) para os municípios mostrados anteriormente na Figura 33:

Município	Arujá	Biritiba-Mirim	Guararema	Itaquaquecetuba	Mogi Das Cruzes	Salesópolis	Suzano
Arujá	0	0	0	1	0	0	1
Biritiba-Mirim	0	0	0	0	1	1	0
Guararema	0	1	0	0	1	0	0
Itaquaquecetuba	1	0	0	0	0	0	1
Mogi Das Cruzes	0	1	0	0	0	0	1
Salesópolis	0	1	1	0	0	0	0
Suzano	0	0	0	1	1	0	0

Vemos que agora nesse caso podemos ter assimetrias: não necessariamente se $c$ é um dos $k$ municípios mais próximos de $l$, $l$ será um dos vizinhos mais próximos de $c$. Isso pode ser visualizado no exemplo apresentado visto que Mogi das Cruzes é vizinha de Guararema (pois está entre suas duas cidades mais próximas) mas Guararema não é vizinha de Mogi, que tem outros dois municípios que são mais próximos (Figura 37).

viz_sp2 <- centroides_coordenadas %>% 
  knearneigh(k=3) %>% 
  knn2nb()

viz_linhas2 <- nb2lines(nb = viz_sp2, 
                       coords = centroides_coordenadas) %>% 
  st_as_sf()
st_crs(viz_linhas2) <- st_crs(malha_sp)

ggplot(malha_sp, aes(geometry=geom)) + 
  geom_sf(fill = "lightgreen") + 
  geom_sf(data = viz_linhas2, aes(geometry=geometry)) + 
  ggtitle("Matriz de vizinhança por k = 3 vizinhos mais próximos, estado de São Paulo.") +
  theme_void()

Figura 37: Matriz de vizinhança por k vizinhos mais próximos dos munícipios do estado de São Paulo.

O script cria uma matriz de vizinhança baseada nos 3 vizinhos mais próximos (k=3) a partir dos centroides dos municípios, desenha as linhas de conexão entre eles e depois plota isso sobre o mapa do estado de São Paulo.

Por distância

Uma terceira estratégia que podemos adotar é a definição de vizinhos a partir de uma distância pré-definida. Nesse caso, a célula $l,c$ será 1 se $c$ estiver até a uma distância $d$ pré-definida de $l$. Para definir essa distância, é importante entender como as distâncias entre as áreas em questão estão distribuídas.

Vamos exibir a matriz de distâncias e sua distribuição para o nosso exemplo de 7 municípios do estado de São Paulo:

	Arujá	Biritiba-Mirim	Guararema	Itaquaquecetuba	Mogi Das Cruzes	Salesópolis	Suzano
Arujá	0,0	40,4	27,1	8,5	24,5	53,4	24,7
Biritiba-Mirim	40,4	0,0	22,2	36,7	17,9	18,9	29,6
Guararema	27,1	22,2	0,0	28,5	20,5	27,9	32,7
Itaquaquecetuba	8,5	36,7	28,5	0,0	19,3	52,0	16,5
Mogi Das Cruzes	24,5	17,9	20,5	19,3	0,0	35,2	13,4
Salesópolis	53,4	18,9	27,9	52,0	35,2	0,0	48,0
Suzano	24,7	29,6	32,7	16,5	13,4	48,0	0,0

Geralmente, trabalhamos com conjunto de áreas muito maiores (estado de São Paulo: 645 municípios), o que torna inviável olhar uma matriz de distância entre cada área (uma matriz 645x645 = 416.025 distâncias). Portanto, é comum recorrermos a métricas do tipo resumo, como médias, medianas e quantis da distribuição de distâncias - ou até mesmo a técnicas gráficas como um histograma - para a definição de um ponto de corte (veremos a seguir). Por enquanto, parece razoável estabelecer o limite de $d=20km$ para definição de nossa matriz de vizinhança com base nas distâncias:

Município	Arujá	Biritiba-Mirim	Guararema	Itaquaquecetuba	Mogi Das Cruzes	Salesópolis	Suzano
Arujá	0	0	0	1	0	0	0
Biritiba-Mirim	0	0	0	0	1	1	0
Guararema	0	0	0	0	0	0	0
Itaquaquecetuba	1	0	0	0	1	0	1
Mogi Das Cruzes	0	1	0	1	0	0	1
Salesópolis	0	1	0	0	0	0	0
Suzano	0	0	0	1	1	0	0

Agora, temos uma matriz simétrica novamente - pois se $l$ está a uma distância de até 20km de $c$, então o contrário também é verdadeiro. Vemos também que o número de vizinhos muda: alguns municípios como Mogi das Cruzes possuem 3 vizinhos pois têm três municípios em até um raio de 20km de distância; outros possuem 1, como Arujá; e Guararema neste caso ficou sem vizinhos, pois está mais afastada dos demais.

Vamos ver como essa situação se aplica ao contexto dos 645 municípios do estado de São Paulo. Primeiro, vamos calcular as distâncias entre os centroides dos municípios. Para isso, convertemos a malha para o sistema de coordenadas UTM (31982), que permite termos coordenadas em metros, ao invés de graus:

# convertendo para UTM, calculando as coordenadas dos centroides novamente e por fim, as distancias:

centroides_coordenadas_metros <- malha_sp %>%
  st_transform(crs = 31982) %>% 
  st_centroid() %>% 
  st_coordinates()

distancias_sp <- centroides_coordenadas_metros %>% 
  dist() 

# um objeto da classe dist
class(distancias_sp)

#> [1] "dist"

#> [1] "dist"
# convertendo para km
distancias_sp <- distancias_sp/1000

distancias_sp %>% 
  as.vector() %>% 
  hist(main = "Histograma das distâncias",
       xlab = "Distância (km)",
       ylab = "Frequência",
       col = "royalblue")

Figura 38: Histograma das distâncias entre os municípios de São Paulo.

O script começa transformando a malha espacial malha_sp para o sistema de coordenadas UTM (EPSG 31982), que utiliza metros como unidade de medida. Em seguida, calcula o centroide de cada polígono (isto é, o ponto central de cada município) e extrai as suas coordenadas X e Y.

distancias_sp <- centroides_coordenadas_metros %>% dist(): A partir dessas coordenadas, é calculada a matriz de distâncias entre todos os centroides, resultando em um objeto do tipo dist, onde as distâncias estão inicialmente em metros.

Depois, todas as distâncias são convertidas de metros para quilômetros, dividindo os valores por 1.000.

Por fim, o código transforma essas distâncias em um vetor e constrói um histograma, visualizando a distribuição da frequência das distâncias (em quilômetros) entre os municípios da malha.

Vemos na Figura 38 que a distribuição das distâncias entre os municípios de São Paulo se encontra na faixa de 200-250km, podendo chegar até 800km. Vamos observar os quantis dessa distribuição:

quantile(distancias_sp,
         probs = c(0.01, 0.025, 0.05, 0.10, 0.25, 0.5))

#>        1%      2.5%        5%       10%       25%       50% 
#>  28.12763  44.34926  64.25050  94.52491 163.10356 260.96230

O comando calcula o valor das distâncias que se encontram nos percentis 1%, 2,5%, 5%, 10%, 25% e 50% da distribuição de distâncias entre os municípios.

É ideal que não escolhamos um valor muito alto - escolher 60km, por exemplo, implicará que 5% de todas as relações entre municípios serão consideradas relações de vizinhança - o que pode ser muita coisa. Vamos tentar com o ponto de corte $d = 40km$:

viz_sp3 <- centroides_coordenadas_metros %>% 
  dnearneigh(d1=0, d2=40000)

viz_sp3

#> Neighbour list object:
#> Number of regions: 645 
#> Number of nonzero links: 8464 
#> Percentage nonzero weights: 2.034493 
#> Average number of links: 13.12248

O script cria uma lista de vizinhança espacial chamada viz_sp3 baseada na distância entre centroides.

• o comando centroides_coordenadas_metros: contém as coordenadas X e Y dos centroides de cada município, em metros (no sistema UTM).

• dnearneigh(d1=0, d2=40000): Cria a vizinhança espacial considerando que d1 = 0 é a distância mínima entre vizinhos é zero metros (ou seja, qualquer distância positiva serve). E d2 = 40000 é a distância máxima para serem considerados vizinhos é 40.000 metros (40 km).

viz_linhas3 <- nb2lines(nb = viz_sp3, 
                       coords = centroides_coordenadas) %>% 
  st_as_sf()
st_crs(viz_linhas3) <- st_crs(malha_sp)

ggplot(malha_sp, aes(geometry=geom)) + 
  geom_sf(fill = "plum1") + 
  geom_sf(data = viz_linhas3, aes(geometry=geometry),
          alpha=0.4) + 
  ggtitle("Matriz de vizinhança por distância (40km), estado de São Paulo.") +
  theme_void()

Figura 39: Matriz de vizinhança por distância dos munícipios do estado de São Paulo.

O script primeiro cria linhas de vizinhança conectando os centroides dos municípios que são considerados vizinhos segundo o objeto viz_sp3. Essas linhas são convertidas para um formato espacial (sf) para que possam ser usadas em visualizações no ggplot2. Em seguida, o sistema de referência espacial dessas linhas (viz_linhas3) é ajustado para ser o mesmo da malha de municípios (malha_sp), garantindo que os dados estejam sobre a mesma base de coordenadas.

Depois, o script monta o mapa: ele desenha os polígonos dos municípios preenchidos com uma cor lilás clara, sobrepõe as linhas de vizinhança com transparência (para dar destaque ao mapa de fundo), adiciona um título explicativo e remove todos os elementos visuais desnecessários (como eixos, grades e bordas) para deixar o mapa mais limpo e focado na informação espacial.

Vemos na Figura 39 todos os municípios possuem vizinhos, mas algumas áreas são bem mais densas do que outras - municípios rurais tendem a ser maiores e mais espaçados entre uns e outros, enquanto os grandes centros parecem possuir municípios menores e mais próximos um dos outros. Isso se confirma quando vemos o número de vizinhos por município (Figura 40):

# extrair numero de vizinhos de cada municipio

num_vizinhos3 <- viz_sp3 %>% 
  lapply(length) %>% 
  unlist()

malha_sp %>% 
  mutate(num_vizinhos = num_vizinhos3) %>% 
  ggplot(aes(geometry=geom)) + 
  geom_sf(aes(fill=num_vizinhos)) + 
  scale_fill_viridis_c(name="Número de vizinhos",
                       breaks = seq(0,30,by=5)) +
  ggtitle("Número de vizinhos de cada município, segundo a estratégia\nde vizinhança por distância (40km)") +
  theme_void()

Figura 40: Número de vizinhos segundo matriz de vizinhança por distância de 40km, munícipios do estado de São Paulo.

O script calcula quantos vizinhos cada município possui com base na lista de vizinhança viz_sp3. Primeiro, ele aplica a função length em cada elemento da lista para contar o número de vizinhos de cada município. Esses valores são então organizados em um vetor simples (num_vizinhos3).

Depois, o script usa a malha espacial dos municípios (malha_sp), adiciona uma nova coluna chamada num_vizinhos contendo o número de vizinhos de cada município, e cria um mapa. No mapa, os municípios são coloridos de acordo com a quantidade de vizinhos, utilizando a escala de cores contínuas viridis, que é perceptível para pessoas com daltonismo e facilita a interpretação. O título explica que o critério de vizinhança foi baseado em uma distância de até 40 km entre municípios, e o tema theme_void() é usado para deixar o mapa limpo, sem eixos nem grades.

Dada as diferenças de resultados obtidos, vemos que a escolha da definição da matriz de vizinhança é importante e pode impactar os resultados obtidos nas análises subsequentes. Por isso, é importante ter conhecimento do agravo em questão e entender quais fatores são determinantes para caracterização de uma relação entre um município e outro (divisão de fronteiras, fluxo migratório e/ou pendular, distância, entre outros). Os tópicos a seguir utilizam as matrizes de vizinhança definidas para aplicação de técnicas de análise espacial aos dados.

Testando a autocorrelação espacial global

Durante o capítulo, comentamos sobre o efeito da correlação espacial - o pressuposto de que o espaço - ou seja, quais municípios estão próximos uns dos outros - influencia na determinação do agravo. Esse pressuposto é importante, visto que determinadas abordagens estatísticas - como modelos de regressão - possuem a independência como pressuposto, e caso assumamos independência entre os dados quando há correlação espacial entre eles podemos obter estimativas enviesadas. Mesmo quando o agravo não é diretamente ou intuitivamente relacionado ao espaço, outros fatores determinantes podem se distribuir espacialmente na população e gerar uma dependência espacial do processo de interesse. Por isso, é sempre importante verificar a distribuição dos dados de maneira exploratória e testar a existência dessa correlação.

Geralmente testamos para a autocorrelação espacial - ou seja, a correlação da variável com ela mesma ao longo do espaço - utilizando índices globais conhecidos, como o Índice I de Moran ou o Índice C de Geary.

# style = "W" corresponde a abordagem informada de somar as arestas e obter 1 
pesos_viz <- viz_sp %>% nb2listw(style = "W") 

moran.test(
  dengue_2020$inc,
  pesos_viz
)

#> 
#>  Moran I test under randomisation
#> 
#> data:  dengue_2020$inc  
#> weights: pesos_viz    
#> 
#> Moran I statistic standard deviate = 18.034, p-value < 2.2e-16
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic       Expectation          Variance 
#>      0.4408150828     -0.0015527950      0.0006017108

geary.test(
  dengue_2020$inc,
  pesos_viz
)

#> 
#>  Geary C test under randomisation
#> 
#> data:  dengue_2020$inc 
#> weights: pesos_viz   
#> 
#> Geary C statistic standard deviate = 13.006, p-value < 2.2e-16
#> alternative hypothesis: Expectation greater than statistic
#> sample estimates:
#> Geary C statistic       Expectation          Variance 
#>       0.544092361       1.000000000       0.001228851

Testando a autocorrelação espacial local

Além de um indicador que estima a existência ou ausência de autocorrelação espacial do processo no conjunto de dados como um todo, pode ser interessante averiguar a existência de fenômenos que ocorrem localmente. Mesmo com ausência de autocorrelação espacial global, podemos observar algumas manifestações locais de autocorrelação, que foram mascaradas ao olhar para a região como um todo. Esses indicadores locais permitem uma exploração mais detalhada das áreas envolvidas em possíveis aglomerados da doença (autocorrelação local positiva), ou a identificação de possíveis outliers (autocorrelação local negativa).

Indicadores Locais de Associação Espacial

O índice de Moran Local (que é um LISA - Local Indicators of Spacial Assocation) é uma adaptação do índice I de Moran para um contexto local, calculado para cada área $i$. Possui interpretação semelhante, mas voltada para cada área, medindo seu grau de associação com suas áreas vizinhas.

lisa_dengue2020 <- localmoran_perm(dengue_2020$inc, pesos_viz) 

head(lisa_dengue2020)

#>            Ii         E.Ii     Var.Ii        Z.Ii Pr(z != E(Ii))
#> 1  1.40116921  0.010362547 0.20960579  3.03784158    0.002382792
#> 2 -0.24463453 -0.010723014 0.04967033 -1.04954983    0.293925134
#> 3  0.47277760  0.006618033 1.05671851  0.45347667    0.650205542
#> 4  0.02506444  0.003231677 0.10332043  0.06792281    0.945847082
#> 5  0.36274715 -0.005694462 0.12166351  1.05630298    0.290829815
#> 6  0.37898181 -0.007048969 0.07344926  1.42438773    0.154334264
#>   Pr(z != E(Ii)) Sim Pr(folded) Sim   Skewness  Kurtosis
#> 1              0.024          0.012  1.0589117 0.5245390
#> 2              0.300          0.150 -1.0392633 1.2729922
#> 3              0.604          0.302  0.8987260 0.7709809
#> 4              0.744          0.372 -1.7630536 3.2617769
#> 5              0.112          0.056 -1.1696210 0.9156516
#> 6              0.008          0.004 -0.8935898 0.6438619

O script calcula o índice de autocorrelação espacial local (LISA) para a incidência de dengue em 2020 usando permutação aleatória.

A função localmoran_perm():

• Recebe dengue_2020$inc (incidência bruta) e pesos_viz (estrutura de pesos espaciais).

• Calcula o índice LISA para cada município, que mede a autocorrelação local — ou seja, identifica clusters de valores altos (hotspots), valores baixos (coldspots) e áreas de transição.

• Utiliza permutações para gerar uma distribuição nula e calcular p-valores para cada município.

O objeto lisa_dengue2020 armazena:

• O valor do índice local de Moran para cada área.

• Estatísticas associadas como p-valor, valores esperados sob aleatoriedade, e variância.

O comando head(lisa_dengue2020) exibe as primeiras linhas, mostrando os resultados dos primeiros municípios analisados.

Agora vemos, para cada município, uma estimativa de I (Ii) e seu p-valor Pr(z != E(Ii)). Podemos representar esse resultado visualmente na Figura 41:

malha_sp_dengue$local_I <- lisa_dengue2020[,1]
malha_sp_dengue$local_I_p_valor <- lisa_dengue2020[,5]

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = local_I)) +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0,
                       name = "I local") +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"))

Figura 41: Índice de Moran Local da incidência de dengue nos municípios de São Paulo, 2020.

O script associa os resultados do LISA (Índice de Moran Local) à malha espacial dos municípios de São Paulo.

Primeiramente:

• A primeira coluna de lisa_dengue2020 (valores do índice local) é atribuída à nova variável local_I em malha_sp_dengue.

• A quinta coluna de lisa_dengue2020 (p-valor associado) é atribuída à nova variável local_I_p_valor.

Em seguida, o script cria um mapa utilizando ggplot2:

• Cada município é preenchido com uma cor baseada no valor do índice local (local_I).

• Utiliza scale_fill_gradient2(), que aplica uma escala de cores: azul para valores negativos, vermelho para valores positivos e branco para valores próximos de zero.

O mapa resultante destaca visualmente áreas com alta ou baixa autocorrelação local na incidência de dengue em 2015.

Vemos que boa parte do estado concentrou valores de I local próximos de zero. Alguns municípios, no entanto, se destacaram por possuírem valores de I bem maiores que zero, o que significa uma forte autocorrelação espacial local.

Uma forma comum de visualizar os resultados de um LISA é através do chamado LISA map. Nele, utilizamos o p-valor obtido para classificar os municípios com autocorrelação espacial local significativa (p < 0,5) ou não significativa (p >= 0,05). Para os significativos, vemos o sinal da associação: se é direta (I local > 0) ou inversa (I local < 0). Por fim, olhamos para a magnitude da variável: se há uma associação espacial local positiva e a taxa no município é alta (acima da média), então há um cluster Alto-alto, pois é uma região onde altos valores costumam se concentrar. Se a taxa é baixa, contudo, trata-se de um cluster Baixo-baixo, pois há uma concentração de municípios com taxas abaixo do esperado. Quando a associação é significativa e negativa, então há os cenários Alto-baixo e Baixo-alto, que ocorre quando o município é um outlier entre seus vizinhos - se ele está com a taxa alta, seus vizinhos tendem a estar baixo e vice-versa.

No resultado do comando, esse agrupamento já é realizado e armazenado no atributo quadr.

quadrantes <- attr(lisa_dengue2020, "quadr")$mean

malha_sp_dengue$quadrante <- case_when(
  quadrantes == "High-High" ~ "Alto-Alto",
  quadrantes == "Low-Low" ~ "Baixo-Baixo",
  quadrantes == "High-Low" ~ "Alto-Baixo",
  quadrantes == "Low-High" ~ "Baixo-Alto"
)

malha_sp_dengue <- malha_sp_dengue %>% 
  mutate(LISA_resultado = case_when(
    local_I_p_valor < 0.05 ~ quadrante,
    local_I_p_valor >= 0.05 ~ "Não significativo"
  ))

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = LISA_resultado)) +
  scale_fill_manual(
    values = c("#e41a1c", "#4daf4a", "#377eb8", "gray90")
  ) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"))

Figura 42: LISA Map da incidência de dengue nos municípios de São Paulo, 2020.

O script organiza os resultados do teste LISA em quadrantes e gera um mapa temático para interpretar os padrões espaciais da incidência de dengue em 2020.

Primeiramente:

• A informação sobre o tipo de associação espacial (quadrante) é extraída dos atributos do objeto lisa_dengue2020, identificando se cada município pertence a “High-High”, “Low-Low”, “High-Low” ou “Low-High”.

• Cria-se uma nova variável quadrante em malha_sp_dengue, traduzindo esses nomes para português.

Em seguida:

• Uma nova variável LISA_resultado é criada:

  • Municípios com p-valor menor que 0,05 são classificados no quadrante correspondente (“Alto-Alto”, “Baixo-Baixo”, etc.).
  • Municípios com p-valor maior ou igual a 0,05 são classificados como “Não significativo”.

O mapa final é construído com ggplot2:

• Cada município é preenchido com uma cor de acordo com o seu resultado no LISA:
  • “Alto-Alto” (vermelho) indica áreas de alta incidência cercadas por outras áreas de alta incidência (hotspots).
  • “Baixo-Baixo” (azul) indica áreas de baixa incidência cercadas por outras de baixa incidência (coldspots).
  • “Alto-Baixo” (verde) e “Baixo-Alto” (verde) indicam áreas de potencial outlier espacial.
  • “Não significativo” (cinza) indica ausência de autocorrelação estatisticamente significativa.
• O layout do mapa é limpo com theme_void(), e a legenda é posicionada abaixo do mapa para melhor organização visual.

Este mapa LISA permite identificar padrões locais de agregação e outliers espaciais na distribuição da dengue.

Com o resultado, temos os agrupamentos dos municípios de acordo com o resultado do LISA map (Figura 42). Vê-se que a maioria dos municípios não apresentou significância estatística; contudo, alguns aglomerados foram identificados - ressalta-se um cluster de baixa incidência no Sul do estado (grupo “Baixo-Baixo”, em azul), e alguns pontos de clusters de alta incidência (grupo “Alto-Alto”, em vermelho), em regiões já identificadas como destoantes nas visualizações anteriores. O LISA map aponta, inclusive, alguns municípios outliers - que estão com valores de baixa incidência, mas rodeados de municípios com taxas mais altas (grupo “Baixo-Alto”, em verde).

Kernel de área

Como vimos no Módulo 2, podemos interpolar uma superfície contínua de densidade a partir de pontos que representam casos de uma doença. Nesse caso, não temos a localização pontual dos casos - apenas a contagem de casos por municípios, ou outra unidade de área. A aplicação da técnica consiste, portanto, na representação dessas áreas como pontos (geralmente se considera o centroide da área ou alguma outra localização de referência, como a sede do município) - e na ponderação desses pontos com a contagem agregada de casos para cada região.

Para essa transformação de dados de área para pontos, vamos utilizar novamente os centroides dos municípios (Figura 43). Além disso, vamos utilizar o pacote patchwork para organizar as imagens de uma forma mais simples. Caso não tenha instalado, siga os comandos abaixo:

# se não estiver instalado, rodar:
install.packages("patchwork")
library(patchwork)

malha_sp_pontos <- malha_sp %>% 
  st_centroid()

g_malha_sp <- malha_sp %>% 
  ggplot(aes(geometry=geom)) + 
  geom_sf(fill = "darkolivegreen1") + 
  ggtitle("Malha dos municípios de São Paulo") +
  theme_void()

g_pontos_sp <- malha_sp_pontos %>% 
  ggplot(aes(geometry=geom)) + 
  geom_sf(color = "darkgreen") + 
  ggtitle("Municípios de São Paulo, transformados\nem pontos (centroides)") +
  theme_void()

g_malha_sp | g_pontos_sp

Figura 43: Divisão territorial e centroides dos municípios do estado de São Paulo.

Agora, temos nossos pontos para aplicação do método Kernel. Utilizaremos o Kernel por atributo - onde cada ponto (cada município) terá um valor que corresponde ao que desejamos visualizar, neste caso, a taxa de incidência de uma doença.

Para isso, precisamos criar um objeto da classe ppp (point pattern).

library(spatstat)
library(nngeo)

## Contorno do estado

contorno_sp <- malha_sp %>% 
  st_transform(crs = 31982) %>% 
  st_union() %>% 
  st_remove_holes()

# Transformando o contorno para janela (owin)
sp_owin <- contorno_sp %>% as.owin()

sp_ppp <- ppp(
  centroides_coordenadas_metros[,1], # Coordenada em X
  centroides_coordenadas_metros[,2], # Coordenada em Y
  sp_owin # Janela
)

plot(sp_ppp, pch = 19, cex = 0.5, main = "Objeto PPP dos municípios")

Figura 44: Visualização do objeto de classe ppp contendo os centroides dos municípios e o contorno do estado de São Paulo.

Agora temos um objeto do tipo ppp - point pattern, tal qual utilizamos no Módulo 2. Agora, vamos visualizar de forma suavizada no espaço os casos de dengue em 2015 no estado de São Paulo. Podemos voltar a utilizar a função density(), para estimação do kernel, mas teremos que nos atentar a dois parâmetros:

• sigma: Largura de banda - indicará a ordem de suavização aplicada aos dados (maior ou menor).
• weights: o parâmetro utilizado para o kernel de atributo. Aqui é importante informarmos o atributo que estamos interessados em suavizar, como por exemplo, a incidência de dengue.
• eps: indicará a resolução (qualidade gráfica) do resultado. Quanto menor, maior resolução.

Vamos testar as seguintes larguras de banda: 5000, 10000, 20000 e 50000.

kernel_den_2015_b5000 <- density(
  sp_ppp,
  5000,
  weights = dengue_2015$inc,
  scalekernel = T,
  diggle = T,
  eps=2000
)

kernel_den_2015_b10000 <- density(
  sp_ppp,
  10000,
  weights = dengue_2015$inc,
  scalekernel = T,
  diggle = T,
  eps=2000
)

kernel_den_2015_b20000 <- density(
  sp_ppp,
  20000,
  weights = dengue_2015$inc,
  scalekernel = T,
  diggle = T,
  eps=2000
)

kernel_den_2015_b50000 <- density(
  sp_ppp,
  50000,
  weights = dengue_2015$inc,
  scalekernel = T,
  diggle = T,
  eps=2000
)

Agora, vamos gerar gráficos para compará-los lado a lado. Podemos transformar o resultado em uma tibble():

kernel_den_2015_b5000 <- kernel_den_2015_b5000 %>% 
  as_tibble()
head(kernel_den_2015_b5000)

#> # A tibble: 6 × 3
#>         x        y    value
#>     <dbl>    <dbl>    <dbl>
#> 1 284117. 7498067. 6.26e-14
#> 2 284117. 7500065. 1.77e-13
#> 3 286114. 7494070. 3.57e-14
#> 4 286114. 7496068. 1.39e-13
#> 5 286114. 7498067. 4.60e-13
#> 6 286114. 7500065. 1.30e-12

Nessa tibble, vemos a coordenada estimada em X, Y e o valor de densidade estimado pelo kernel. Lembramos que o kernel sempre estima a densidade do evento em cada localização.

g_kernel_5000 <- kernel_den_2015_b5000 %>% 
  ggplot() + 
  geom_tile(aes(x=x, y=y, fill=value)) + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "sigma: 5000", fill = "densidade") +
  coord_fixed() + theme_void()

kernel_den_2015_b10000 <- kernel_den_2015_b10000 %>% 
  as_tibble()

g_kernel_10000 <- kernel_den_2015_b10000 %>% 
  ggplot() + 
  geom_tile(aes(x=x, y=y, fill=value)) + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "sigma: 10000", fill = "densidade") +
  coord_fixed() + theme_void()

kernel_den_2015_b20000 <- kernel_den_2015_b20000 %>% 
  as_tibble()

g_kernel_20000 <- kernel_den_2015_b20000 %>% 
  ggplot() + 
  geom_tile(aes(x=x, y=y, fill=value)) + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "sigma: 20000", fill = "densidade") +
  coord_fixed() + theme_void()

kernel_den_2015_b50000 <- kernel_den_2015_b50000 %>% 
  as_tibble()

g_kernel_50000 <- kernel_den_2015_b50000 %>% 
  ggplot() + 
  geom_tile(aes(x=x, y=y, fill=value)) + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "sigma: 50000", fill = "densidade") +
  coord_fixed() + theme_void()

(g_kernel_5000 | g_kernel_10000) /
  (g_kernel_20000 | g_kernel_50000)

Figura 45: Estimativa de densidade de kernel da incidência de dengue por município no estado de São Paulo, 2015.

Vemos na Figura 45 que a utilização de um sigma baixo, como $sigma=5000$, gerou uma visualização muito concentrada ainda no centroide dos municípios, suavizando muito pouco e dificultando a percepção de padrões espaciais. Com valores de largura de banda maiores, como $sigma=10000$ ou $sigma=20000$ vemos uma suavização interessante da incidência, permitindo ver regiões de maior concentração da doença naquele ano. Vemos o destaque principalmente na região Noroeste de São Paulo, marcada em tons vermelho na escala de cores. O $sigma=50000$ já realizou uma suavização demasiadamente grosseira, permitindo ver a distribuição em termos macro apenas.

Perceba também, que a escala de densidade muda em cada gráfico: com larguras de banda menores, os tons vermelhos representam uma densidade muito mais alta (9e10-6 = $0,000009$), e na maior largura de banda os tons vermelhos representam uma densidade menor (1.5e10-6 = $0,0000015$, seis vezes menor).

Vamos comparar a visualização que tínhamos anteriormente da incidência bruta com a visualização suavizada por kernel na Figura 46:

# gráfico feito anteriormente
g_tradicional <- ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc)) +
  labs(title = "Incidência bruta", fill = "Incidência (por 10.000 hab.)") + 
  scale_fill_distiller(palette = "Reds", direction=1) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"))

g_kernel <- malha_sp %>% 
  st_transform(crs = 31982) %>% 
  ggplot() + 
  geom_tile(
    data = kernel_den_2015_b10000,
    aes(x=x, y=y, fill = value)
  ) +
  geom_sf(aes(geometry=geom), fill="transparent") + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "Estimativa de densidade de kernel, Sigma=10000", fill = "Densidade") +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"))

g_tradicional | g_kernel

Figura 46: Comparação da visualização dos valores burtos e da estimativa de densidade de kernel da incidência de dengue nos municípios do estado de São Paulo, 2015.

Vemos que, por meio do Kernel por atributo, conseguimos obter uma superfície contínua de densidade de incidência, sobre a qual podemos capturar possíveis padrões espaciais existentes nos dados. Na Figura 46, ao olhar a distribuição da densidade de incidência de dengue em 2015 no estado de São Paulo, notamos a concentração de valores no tom vermelho mais escuro (e mais alto) da escala no Noroeste do estado, além de outros dois ou três aglomerados que atingem os tons amarelados e chamam atenção. Contudo, ao olharmos um dado suavizado por Kernel, olhamos para densidade - o que dificulta a interpretação dos resultados obtidos em termos de casos ou incidência. Veremos a seguir um método de suavização que mantém a mesma unidade da variável de interesse, ou seja, conseguimos suavizar os dados e interpretá-los em termos de incidência.

Método Bayesiano Empírico

Ainda com o intuito de suavização de taxas no mapeamento de doenças, o Método Bayesiano Empírico introduz uma média das taxas que é ponderada de acordo com o tamanho da população do município. Essa média pode ser em relação ao valor global (taxa no estado como um todo, por exemplo) ou local (taxa média nos municípios vizinhos). O método local costuma ser mais interessante no mapeamento de doenças ao fornecer uma suavização mais realista - com base no contexto específico do município e seus vizinhos ao invés de uma média generalizada. A fórmula do método é a seguinte:

Perceba que se trata de uma operação de média ponderada entre dois valores ($r_i$ e $\mu_i$), onde o peso é dado por $w_i$.

• $\hat{\theta}_i$ é a taxa suavizada no município $i$.
• $r_i$ é a taxa bruta de incidência no município $i$.
• $\mu_i$ é a taxa média nos vizinhos de $i$ (ou taxa média geral, se for utilizado o método Bayesiano Empírico Global).
• e $w_i$ é o peso a ser considerado nessa operação de média.

O peso $w_i$ é calculado com base nas variâncias da taxa global e da taxa no município $i$. Quanto maior a população, menor é a variância e assim $w_i$ é próximo de 1, dando mais peso à taxa bruta. Quanto menor a população, temos uma variância maior e taxas mais instáveis, e assim mais peso é dado para a média dos vizinhos.

No R temos uma função que realiza esses cálculos diretamente - EBLocal(), do pacote spdep. Vamos aplicá-la ao nosso mesmo caso de incidência de dengue em São Paulo em 2015:

# Relembrando a estrutura dos dados
head(dengue_2015)

#>   cod_ibge cod_mun6               nome_mun  ano   pop casos       inc
#> 1  3500105   350010             Adamantina 2015 34285  1140 332.50693
#> 2  3500204   350020                 Adolfo 2015  3903   121 310.01793
#> 3  3500303   350030                  Aguaí 2015 32192  2487 772.55219
#> 4  3500402   350040         Águas da Prata 2015  7558   336 444.56205
#> 5  3500501   350050       Águas de Lindóia 2015 17690   165  93.27304
#> 6  3500550   350055 Águas de Santa Bárbara 2015  6313    16  25.34453

A função EBLocal() recebe três parâmetros: o número de casos em cada área, a população de cada área, e a estrutura de vizinhança, que definimos anteriormente. Ou seja, vamos optar por uma das estratégias experimentadas (vizinhança por conectividade, vizinhos mais próximos, ou por distância).

inc_bayes_empirico <- EBlocal(
  ri = dengue_2015$casos, 
  ni = dengue_2015$pop,
  nb = viz_sp # matriz de vizinhança definida por contiguidade
)

head(inc_bayes_empirico)

#>           raw         est
#> 1 0.033250693 0.033242118
#> 2 0.031001793 0.030968000
#> 3 0.077255219 0.077160864
#> 4 0.044456205 0.044903289
#> 5 0.009327304 0.009712328
#> 6 0.002534453 0.002692929

Vemos que o objeto retorna as taxas bruta (raw) e suavizada (est). Vamos extrair a segunda coluna (suavizada) e multiplicá-la por 10.000, como fizemos com a incidência:

dengue_2015 <- dengue_2015 %>% 
  mutate(inc_EB = 1e4*inc_bayes_empirico[,2])

malha_sp_dengue$inc_EB <- dengue_2015$inc_EB

malha_sp_dengue %>% 
  pivot_longer(cols = c(inc, inc_EB), names_to = "metodo", values_to = "inc") %>%
  mutate(metodo = factor(metodo, 
                         levels = c("inc", "inc_EB"),
                         labels = c("Incidência bruta", "Incidência suavizada por EB Local"))) %>% 
  ggplot(aes(geometry=geom)) + 
  geom_sf(aes(fill = inc)) +
  labs(fill = "Incidência (por 10.000 hab.)") + 
  scale_fill_distiller(palette = "Reds", direction=1) +
  facet_wrap(~metodo) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"),
        plot.title = element_text(hjust=0.5, size=16))

Figura 47: Incidência bruta e suavizada de dengue nos municípios de São Paulo, 2015.

Vemos que outros padrões chamam a atenção ao comparar a taxa bruta e a taxa suavizada na Figura 47: vê-se que houve redução de algumas flutuações principalmente no sul do estado, e alguns municípios tiveram uma elevação na taxa após a suavização. As maiores diferenças entre incidência bruta e suavizada se dão principalmente nos municípios menores - onde a baixa população implica num menor peso para a incidência bruta calculada, levando a uma maior influência das taxas dos vizinhos.